Formation and evolution of clumpy tidal tails around globular 

clusters 

R. Capuzzo Dolcetta 
dolcetta@uniromal . it 
^ ' P. Di Matteo 

o i 

^i p.dimatteo@uniromal.it 

and 

P. Miocchi 

miocchi@uniromal . it 

Dep. of Physics, Universitd di Roma La Sapienza, P.le Aldo Moro 2, 00185 Rome, Italy 



ABSTRACT 



> 

CO 

CO 

\o 
o 
"*t 

o 

We present some results of numerical simulations of a globular cluster orbiting 
in the central region of a triaxial galaxy on a set of 'loop' orbits. Tails start 
-O ■ forming after about a quarter of the globular cluster orbital period and develop, 

in most cases, along the cluster orbit, showing clumpy substructures as observed, 
for example, in Palomar 5. If completely detectable, clumps can contain about 
7OOOM each, i.e. about 10% of the cluster mass at that epoch. The morphology 
of tails and clumps and the kinematical properties of stars in the tails are studied 
and compared with available observational data. Our finding is that the stellar 
velocity dispersion tends to level off at large radii, in agreement to that found 
for Ml 5 and u Centauri. 

Subject headings: methods: n-body simulations, globular clusters: general, galax- 
ies: kinematics and dynamics 



1. Introduction 

Since Shapley's pioneering work (Shapley 1918), globular clusters (GCs) have played 
a key-role in our understanding of the Universe and of the manner in which our Galaxy 
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formed: in the Milky Way they are the oldest stellar systems found, with ages in the range 
12 to 15 Gyr, so to represent tracers of the early formation history of the Galaxy. 
They are the best systems to study stellar dynamics, having relaxation times smaller than 
their age, so that, at least in the core, stars are expected to have lost memory of their initial 
conditions (Binney and Tremaine 1987). In the early 1980s a number of approximated 
numerical studies of spherical self-gravitating systems (Cohn 1980) showed that the central 
density tends to increase dramatically over the time, so that ultimately a central power-law 
cusp is produced in the central region, even if these systems have an early evolutionary 
phase that resembles the King sequence of cluster models (King 1966). Together with the 
slow collapse of the core, star evaporation occurs: the approach to equipartition implies that 
the more massive stars sink toward the center of the cluster, while lighter stars expand their 
orbit. Core collapse can be halted by the presence of hard binaries which, acting as energy 
sources, heat the central core by 3-bodies encounters (Henon 1961; Ostriker 1985). 
Internal processes are not the only responsible of dynamical evolution in these systems: 
perturbations due to an external field (in particular, shocks due to the passage through the 
galactic disk and to the interaction with the bulge) can accelerate significantly the evolution 
of a globular cluster. Indeed, it is commonly accepted that the present globular cluster 
population represents the survivor of an initially more numerous one, depauperated by many 
disruptive processes (Murali & Weinberg 1997a,b; Fall & Zhang 2001). 
There is observational evidence that the globular cluster system (GCS) radial profile is less 
peaked than that of halo stars in our galaxy, M31 (Capuzzo Dolcetta & Vignola 1997), 
M87 and M49 (Grillmair et al. 1986; McLaughlin 1995), as well as in three galaxies of 
the Fornax cluster (Capuzzo Dolcetta & Donnarumma 2001) and of 11 elliptical galaxies 
(Capuzzo Dolcetta & Tesseri 1999). This fact leads to the hypothesis (Capuzzo Dolcetta & 
Tesseri 1997) that the two systems (halo and GCS) originally had the same profile and that, 
afterwards, the GCS evolved mainly due to two complementary effects: tidal interaction 
with the galactic field (which causes less concentrated clusters to disintegrate more rapidly) 
and dynamical friction (which induces massive globular clusters to decay in the central 
galactic region in less than 10 8 years, see Capuzzo Dolcetta & Vicari (2003)). External 
tidal fields have the effect of inducing the evolution of the shape of the mass function of 
individual clusters, because of the preferential depletion of low-mass stars (Baumgardt & 
Makino 2003) as a consequence of two-body relaxation . Strong evidence that the tidal field 
plays a fundamental role in the evolution of mass functions was achieved by the discovery 
that their slopes correlate more strongly with the cluster location in the Milky Way than 
with the cluster metallicity (Djorgovski et al. 1993). 

In the last decade, many observational evidences of the interaction of GCs with the tidal 
field have been found. Firstly, Grillmair et al. (1995), using colour- magnitude selected star 
counts in a dozen of galactic GCs, showed that in the outer parts of these clusters the stellar 
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surface density profiles exceeded the prediction of King models, extending also outside the 
tidal radius of the corresponding King model. Other results confirmed Grillmair's findings 
(Lehmann & Scholz 1997; Testa et al. 2000; Leon et al. 2000; Siegel et al. 2001; Lee et al. 
2003); all these works suggest that many GCs are likely surrounded by haloes or tails, made 
up of stars which were tidally stripped from the system. This was the state of the art until 
the spectacular findings of two tidal tails emanating from the outer part of the Palomar 5 
globular cluster and covering an arc of 10 degrees on the sky, corresponding to a projected 
lenght of 4 kpc at the distance of the cluster (Odenkirchen et al. 2001, 2003), obtained in 
the framework of the Sloan Digital Sky Survey Project (see also http://www.sdss.org). 
One of the relevant observational features of Palomar 5 is the presence of well defined clumps 
in the star distribution along the tails (Odenkirchen et al. 2001, 2003). Also NGC 6254 and 
Palomar 12 seem to show clumpy structures in their tails (Leon et al. 2000). This still 
deserves an exaustive interpretation. Actually, the simulations of Combes et al. (1999) show 
the presence of small clumps (containing about 0.5% of the total number of stars of the 
cluster) in the tidal tails. The authors attribute the formation of these clumps to strong 
gravitational shocks suffered by the cluster. On another hand, Dehnen et al. (2004) were 
not able to reproduce the clumps observed in the tails of Pal 5, even adopting a realistic 
galactic model, so that they argued that these structures could be due to the effect of 
Galactic sub-structures not accounted in their simulations (giant molecular clouds, spiral 
arms, dark-matter sub-halos or massive compact halo objects). 

With the aim of understanding better the mechanism of interaction of GCs with the 
external field, in particular with the bulge, and to investigate on the presence of clumps 
in tidal tails, we performed numerical simulations of globular clusters in orbit in a triaxial 
galaxy, aiming also at clarifying the morphological connection between the clusters tidal tails 
and their orbits. 

In the next sections we show the results for globular clusters on 'loop' orbits in an inner region 
of a triaxial galaxy. In particular, in Sect. 2, an overview of the numerical methods adopted 
to perform the simulation are discussed; in Sect. 3 the galaxy and cluster model adopted are 
presented; in Sect. 4 and Sect. 5 we deal with the main results of our work, especially that 
concerning the formation of tidal tails around the cluster and their orientation respect to the 
cluster orbit (Sect.4.1), the radial density profiles of the cluster, as they evolve with time, 
and the presence of clumpy regions in the tails (Sect. 4. 2), the velocity dispersion of stars 
in the cluster (Sect. 4. 3), the estimate of the mass loss rate (Sect. 5.1) and the evolution of 
the global mass functions (Sect. 5. 2); in the last section all the findings are summarized and 
discussed. 
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Numerical method 



All the simulations were performed by means of an implementation of a tree- code car- 
ried out mainly by one of the authors (P.M.). It is based on the algorithm described in 
Barnes & Hut (1986) and adopts multipolar expansions of the potential truncated at the 
quadrupole moment. It was parallelized to run on high performance computers via MPI rou- 
tines, employing an original parallelization approach (Miocchi & Capuzzo-Dolcetta 2002). 
The time-integration of the 'particles' trajectories is performed by a 2 nd order leap-frog algo- 
rithm. This latter uses individual and variable time-steps according to the block-time scheme 
(Aarseth 1985; Hernquist & Katz 1989), in addition with corrections implemented in order to 
keep the same order of approximation also during the time-step change. The maximum time- 
step allowed is At max = 0.01t c (where t c ~ (r^ ore /GM) 1//2 is the core-crossing time of the GC, 
being M its mass and r core the core radius), while the minimum is At min = At max /2 8 , thus 
fastest particles may have a time-step as small as ~ 4 x 10~ 5 £ c . The best criterion we found 
for choosing the time-step of the i-th particle is via the formula: Hiinj^j/aj) 1 / 2 , di/vi}/ 20, 
where Vi is the velocity of the particle relative to its first neighbour, di the distance from its 
first neighbour and o» the modulus of the acceleration. 

To avoid instability in the time-integration, we smoothed the l/r^ interaction potential 
by substiuting l/r„ with a continuous /9-spline function that gives an exactly Newtonian 
potential for r^ > e and a force that vanishes for r^ — > (Hernquist & Katz 1989). In 
all the runs we set e = 1.4 x 10~ 3 r core , so to have (e 3 /GM) 1 / 2 ~ At m i n . Note that such 
value of e is much less than the typical interparticle distance. As regards the quality of the 
orbits time-integration, we checked that the upper bound of the relative error in the energy 
conservation (AE/E) is 10~ 8 per time-step, even in absence of the external field. 



3. Cluster and galaxy models 

3.1. Galaxy model 

The external galactic field due to the bulge is represented by the potential of the 
Schwarzschild model (Schwarzschild 1979). The Schwarzschild model is a non- rotating, self- 
consistent triaxial ellipsoid with axis ratios 2 : 1.25 : 1, typical of many galaxies (Bertola et 
al. 1991). 
Defining adimensional units as 



x' = x/r b , y' = y/r b , z' = z/r b , (1) 
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Tb being the bulge core radius, the potential $(x', y', z') is expressed as the sum of a spheri- 
cally simmetric term, $ r /(r'), (r' = r/rb), which corresponds to the potential given by a den- 
sity distribution following the modified Hubble's law p{r') — p [l + (r') ] , (po = Mb/rf), 
plus two spherical armonics, $i(V,r') and ^2{x',y',r'): 



$(x', y', *') = A [$ r , (/) + $! (z', r') + $ 2 (x', y', r')] , 



where 
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and Mb is the bulge mass. 

The coefficients a have the values: a = 0.06408, c 2 = 0.65456, c 3 = 0.01533, c 4 = 0.48067 

(de Zeeuw & Merritt 1983; Pesce et al. 1992). They have been determined so to have density 

axial ratios roughly constant with r. 

Following Pesce et al. (1992), we will consider Mb = 3 x 10 9 M and r^ = 200 pc, but the 

results obtained in adimensional variables (see Appendix A) are scalable in terms of r^ and 

Mb, for given initial conditions for positions, velocities and rrii/Mb. 



3.2. Cluster model 

As initial cluster model, we chose a multimass King distribution (King 1966; Da Costa 
& Freeman 1976), with 10 mass classes, ranging between 0.12 and 1.2M Q and equally spaced 
in a logarithmic scale. To find the distribution function for each mass class, we integrated 
the Poisson's equation as described in the Appendix B. 

The initial cluster mass function is chosen in the Salpeter form (Salpeter 1955), i.e. dN/dm oc 
m~ 2 ' 35 (see Appendix C for a discussion about the remnants of progenitor stars more mas- 
sive than 1.2M ). We included mass segregation in the initial conditions of our cluster 
model because we wanted to simulate a dynamically-relaxed cluster (supposed to be suf- 
ficiently massive to have frictionally decayed in the central galactic region). The "initial" 
mass of the system is M tat = 3 x 10 5 M = 10~ A Mb, the initial concentration parameter is 
c = log(r t /r corc ) = 1.1 and the central velocity dispersion is a = 9.4 km s _1 = 0.036 xrb/t cross , 
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being t cross = (r^/GM^,) 1 / 2 the bulge crossing time. 

Then, the system was represented by a number (N) of 'particles' lower than the number 

of stars in the real cluster, with masses properly rescaled such to give a total mass equal 

to M tot . The cluster moves on the y-z coordinate plane, following loop orbits of different 

ellipticity 

Ra + Rp 

being R a and R p , respectively, the apocenter and pericenter distances (see Table 1 for orbital 
parameters) . 



4. Results on tidal tails 

In the following, we present the main findings of our work concerning with tidal tails 
structure and evolution. When we refer to the center-of-density of the cluster, we mean a 
mass density- weighted center as defined by Casertano & Hut (1985). 



4.1. Tidal tails formation and morphology 

In all the simulations performed, the cluster starts moving around the galaxy center in 
a clockwise direction (seen from the positive x axis). The different loop orbits have been 
followed for about 30 t croS s- 

In Fig.l, 2, 3, 4, 5, 6, 7 and 8, the formation and subsequent development of tails around 
the globular cluster is shown. 

After about 8 t cross , tidal tails are clearly formed. They continuosly accrete by stars leaving 
the cluster, so that after 30 t croS s, in the case of quasi-circular orbit (e = 0.03), they are 
elongated for more than 3 r& each and contain about 75% of the initial cluster mass. As it 
is clearly visible from these figures, the degree of elongation of the tails along the cluster 
orbital path strongly depends on the ellipticity e (Eq.7) of the cluster orbit. Indeed, while 
in the case of the quasi-circular orbit tails are a clear tracer of the cluster path, in the most 
eccentric orbit (e = 0.57), tails are strictly elongated along the orbital path only when the 
cluster is near the perigalacticon, while at the apogalacticon they tend to deviate from the 
cluster path. Nevertheless, in Miocchi et al. (2004) a remarkable tails — orbit alignment is 
found for clusters moving on quasi-radial orbits in the same bulge potential. However, it is 
important to stress that, in order to perform accurate predictions of the cluster orbit from 
observational detections of tidal tails, it is necessary to look at the spatial distribution of 
stars well outside the cluster (typically 2 — 3 times the cluster limiting radius). Indeed, in 
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the vicinity of the globular cluster, stars in the tails are not aligned with the cluster orbit, 
neither in the case of small ellipticity (see Fig. 9), but they distribute along the peculiar 
S— shape profile not aligned along the orbit. 

The orbit ellipticity also influences the similarity between the two cluster tails. For the 
quasi-circular orbit, these structures are simmetric for the whole duration of the simulation, 
being elongated, at a given time, for the same length. For more eccentric orbits, the leading 
tail tends to be more elongated than the trailing one when going from the apocenter to the 
orbital pericenter and, viceversa, it is less elongated than the trailing tail when the cluster 
moves towards the apocenter. In any case, the tail that precedes the cluster extends always 
slightly below the orbit while the trailing one lies slighlty above this latter, in agreement 
with what observed for Palomar 5. 

The shape and orientation of the tails can be easily understood in the case of a cluster 
moving on a circular orbit in an axysimmetric external field, using a rotating frame of ref- 
erence with the origin in the baricentre of the cluster, with the X-axis pointing towards the 
galactic center, the F-axis parallel to the direction of motion of the cluster and the Z-axis 
orthogonal to the orbital plane. In this reference frame, the galactic tidal field tends to 
accelerate stars along the ±X directions (Heggie & Hut 2003), making stars to escape from 
the system through the Lagrangian points L\ and L2 (which are the two equilibrium points 
located along the X-axis). But the Coriolis acceleration tends to align escaping stars along 
the direction of motion of the cluster around the galaxy, this yelding the peculiar .S-shape 
just outside the cluster, in the inner part of the tails. 



4.2. Density profiles 

In order to describe the tidal debris and to compare our findings with observations 
(Lehmann & Scholz 1997; Testa et al. 2000; Leon et al. 2000; Siegel et al. 2001; Lee et al. 2003; 
Odenkirchen et al. 2003), we studied the radial profile of the volume and surface densities 
(azimutally averaged) as a function of the distance from the cluster center. Obviously, this 
description does not take into account the fact that stars lost from the cluster are not placed 
in a spherically symmetric structure, but it has the advantage to provide a global study of 
both the cluster and the tails that can be easily compared with observational data. In Fig. 10 
and 11, the volume density of the system is shown at different epochs, for the various orbits. 
Once the tails have completely developed, outside the S'-shape distribution, density clumps 
appears. They are symmetrically located in the two tails, as shown in Fig. 12 for the cluster 
on quasi-circular orbit: in this case, the most prominent clumps are located at a distance 
from the cluster center between 0.25r{, and 0.4r^. The density profiles are very similar to 



that of Palomar 5, where clumps are visible in the outer part of the cluster (Odenkirchen et 
al. 2003). 

Of course, the possibility to detect observationally these clumps is strongly related to the 
cluster position along its orbit. Indeed, we computed the contrast density ratio p c i/p*, where 
Pd is the local volume density in the clumps and p* is the background (i.e. the bulge) density 
around them. This ratio is maximum when the cluster is near apogalacticon and decreases 
when moving towards perigalacticon, as shown in Table 2. This is due to two complementary 
effects: when the cluster is near apogalacticon, p* is minimum (according to the galaxy model 
described in Sect. 3.1) and, at the same time, the elongation of the tails tends to compress 
respect to that at perigalacticon (see, for example, Fig. 7 and Fig. 8) and so p c \ increases. If 
completely detectable, clumps can contain about 7000M Q each (i.e. about 10% of the cluster 
mass at that epoch), as in the case of the cluster moving on the quasi-circular orbit after 
SO/ 

In order to study the mass distribution along the tails, we have also evaluated the "linear" 
density for the whole system in the quasi-circular orbit around the galaxy. This study is 
particularly well fitted to investigate the mass distribution because tails form a long and 
thin structure. The upper panel of Fig. 13 shows the linear mass density as a function of the 
curvilinear abscissa s along the system: the absolute maximum in the plot corresponds to 
the cluster location, while the two simmetric relative maxima correspond to the two main 
clumps. These clumps result to be unbound structures (see also Di Matteo et al. (2004)); we 
followed the motion of stars that at a certain time stay in the two clumps: they crowd in the 
clumps for some time and then move away in the outer parts of the tails. Once moved away 
from clumps, these stars tend to disperse along the cluster tails. Also the simmetric location 
of these two clumps respect to the cluster center makes improbable that these structures 
can be due to local disomogenities in the gravitational field along the tails. More probably, 
these clumps are related to cinematical properties of stars in their surroundings. The bottom 
panel in Fig. 13 shows the derivative of the stellar tangential velocity component with respect 
to the curvilinear abscissa s defined above. As is evident, the two clumps correspond to a 
region where this derivative has a negative minimum, which is also the global minimum over 
the whole extension of the tails. This implies that the local velocity of the stars decreases as 
they are approaching the clumps, thus leading to the local overdensity which originates such 
structures. However, the mechanism at the basis of the formation of these structures still 
requires further and more detailed investigationsd that we postpone to next papers. See, 
however, the discussion in Miocchi et al. (2004) for the case of clusters in quasi-radial orbits. 
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4.3. Velocity dispersion 

For the two galactic globular clusters M15 and u Centauri there are observational evi- 
dences that the stellar velocity dispersion remains constant at large radii (Scarpa et al. 2003; 
Drukier et al. 1998). 

Three hypotheses have been raised to justify these findings: 1) tidal heating, as suggested 
by Drukier for M15 (Drukier et al. 1998); 2) the presence of a dark matter halo surrounding 
the clusters (Carraro & Lia 2000); 3) a breakdown of Newton's law of gravity in the weak 
acceleration regime (Scarpa et al. 2003). 

We studied the velocity dispersion profile of our simulated cluster as it would be detected 
if the system was seen along a line-of-sight perpendicular to the cluster orbital plane. In 
Fig. 14, line-of-sight velocities of members of the cluster are plotted versus distance from the 
center, at four different epochs. At t — 0, the velocities decrease moving from the center of 
the cluster outwards, as it is expected from a King model with a tidal cutoff. As the system 
moves through the galaxy and loses stars, the velocity profile varies significantly: it tends to 
decrease until a limiting value and then increases again. This behaviour is very similar to 
that found in M15 (cfr Fig. 8 in Drukier et al. 1998). 

In Fig. 15 the line-of-sight velocity dispersion profile is shown. It is evident from the Figure 
that in the outer part of the cluster, the dispersion tends to level off. This region corresponds 
to that characterized by a power-law volume density profile (see Fig. 10 and 11). Stars in 
this region are escaping from the system and their motion is mostly oriented along the radial 
direction towards the galaxy center. Once escaped, they move around the galaxy weakly 
interacting with each other, with similar orbital parameters, so that the velocity dispersion 
found is coherent with that of a set of particles moving in the triaxial potential adopted. 
The second relevant finding is the decreasing of the velocity dispersion in the inner part 
of the cluster, which could be explained by the quick revirialization of the inner part as 
stellar mass is being lost. This in accordance with the very low velocity dispersion of Pal5 
(Odenkirchen et al. 2002), a cluster which has suffered a great mass loss, as it is now well 
estabilished. 



5. Results on mass loss 

5.1. Mass loss 

To estimate the mass loss from the cluster, we decided to use an 'observational' defini- 
tion. At any given time we compare the cluster local density p gc with the background stellar 
density p*, assuming that a star is actually belonging to the cluster if it is located in a region 
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dense enough to make it distinguishable from the background, i.e. if 

^ > 1, (8) 

P* 

being 

Ap = (p* + Pgc) ~ P*- (9) 

The limiting radius r^ is then defined as the radius of the sphere (centered in the cluster 
density center) in which the cluster 'emerges' from the stellar background. In Fig. 16, the 
evolution of the cluster mass, expressed in units of the initial mass Mo, is shown versus time 
for all the four simulations performed. 

In the case of a cluster moving on orbits with apocenter < 3.5rj,, the mass loss is dramatic: 
after about 30 t croS s the cluster loses about 75% of its mass; the best fit of the mass evolution 
as a function of time is given by: 

^|y = 0.77e"*/ 12 + 0.21, (10) 

where t is expressed in units of the bulge crossing time t cross . In the remaining case, when 
the orbit extends up to 7.5 r^, the mass loss rate considerably diminishes and the cluster 
mass, after 30 t cr0 ss, is still about 60% of its initial value. As is evident in Fig. 16, in this case 
the mass loss rate increases every time the cluster passes at the minimum distance from the 
galaxy center and not all particles which become unbound at perigalacticon are still so while 
moving again to apogalacticon. It is possible to point out a region around the galaxy center 
inside which the cluster suffers more of mass loss: in our case (cluster concentration equal 
to 1.1) this region corresponds roughly to r <4 r&. This conclusion is accordance to what 
found in Miocchi et al. (2004), where great mass loss occurs for clusters with comparable 
central density moving on quasi-radial orbits within such region. 

Finally, we want to stress that the choice of performing some of the simulations with a 
reduced number of particles (N = 1.6 x 10 4 ) did not affect the mass loss rate over the time 
interval of 30 t CT oss- Actually, Fig. 16 clearly shows that, for the cluster in a quasi-circular 
orbit, the mass loss rate is the same using either N = 1.6 x 10 5 (solid line) or a ten times 
smaller N (dashed line). 



5.2. Mass segregation and mass function 

As explained in Sect. 3. 2, we adopted mass segregation in the globular cluster initial 
conditions, for we aim at simulating a dinamically evolved cluster, whose orbit had decayed 
in the inner galactic region due to dynamical friction. As the cluster begins to lose stars, the 
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distribution of stars of different masses in the system starts evolving. This is shown in the 
left column of Fig. 17, where the mean mass of stars populating three different spatial regions 
versus time is plotted. The first region corresponds to the sphere with radius r = 0.016r^ 
(corresponding to r = 3.22 pc, with our choice of Mf, and rj,) centered on the cluster, which 
initially contains 40% of the total mass of the system; the second region is the spherical 
shell with inner and outer radius r = 0.016rt, and r = 0.036rb (r = 7.18 pc), which initially 
contains 80% of the cluster mass; the third region is that outside r = 0.0367V As time passes 
by, low mass stars begin to escape from the system, and the mean stellar mass in the two 
inner regions starts to rise, while the external one remains quite constant. The increasing 
of the mean stellar mass versus time in the central cluster regions is particularly evident in 
the case of the quasi-circular orbit and of the loop with ellipticity e = 0.27, because of the 
greater mass loss in these two cases. Plotting the mean stellar mass in the three regions 
defined above as a function of the fraction of mass lost, we see (Fig. 17, right column) that 
the evolution of the mean mass depends mostly on the fraction of mass lost from the system 
rather than on the number of stars populating the cluster and on the cluster orbital path. 
The differential mass loss obviously influences the shape of the mass function at different 
times. In Fig. 18 the mass function of stars belonging to the cluster is shown at three different 
epochs, when the cluster has lost respectively the 20%, the 35% and the 75% of its initial 
mass. As the cluster loses stars in the galactic field, the mass function evolves towards 
flatter configurations, because of the preferential loss of low- mass stars, that, accordingly to 
the initial mass segregation, are located mostly in the external regions of the cluster. The 
evolution of the mass function appears to be driven by the fraction of mass loss, rather than 
by other parameters (like the total number of stars in the system and the orbital type of 
the parent cluster) confirming the findings of Baumgardt & Makino (2003). This is evident 
from the fact that, for a given fraction of mass lost, the curves found for the different orbits 
in practice coincide. 



6. Conclusions 

The main results of our work may be resumed as follows: 

. Stars are lost from the system along a direction which results from the composition of 
the direction towards the galactic center and the cluster velocity around the galaxy, 
thus leading to the peculiar S-shape found in the outermost region of the cluster. Once 
formed, tidal tails are elongated such to remain parallel to the cluster orbit, with a 
trailing tail that lies slightly inside the orbit and a leading tail slightly outside it. Tails 
are excellent tracers of the cluster orbit near the pericenter, while, at the apocenter, 
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they tend to deviate from the orbital path. 

2. Tidal tails have a clumpy structure which cannot be associated with an episodic mass 
loss or tidal shocks with galactic compact sub-structures, since stars are lost from the 
cluster continuously and since the interaction with the bulge is not episodic. These 
clumps are not bound self-gravitating systems, they are rather due to a local deceler- 
ation of the motion of the stars along the tails. 

3. The observational evidence found for M15 and uo Centauri that the velocity disper- 
sion increases and then remains constant at large radii is explained in terms of the 
so-called 'tidal heating': the stars that evaporate outside the tidal radius of the clus- 
ter undergo mainly the interaction with the external field, thus acquiring the higher 
velocity dispersion pertaining to that field. 
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A. Adimensionalization of the equations 

The equations of motion of the j—th star of the cluster, interacting with all the other 
cluster members and with the bulge are: 

^ = Yl -^ ( r * ~ r ^ + VUb v*iw) ( A1 ) 

where Yj = (xj,yj,Zj) is the position vector, m is the distance between the i— th and 
the j—th particle and U b is the bulge potential. In the case of the Schwarzschild model 
(Schwarzschild 1979): 



U b (x,y,z) = AnGM, 
-3c 3 



1, (r r- r^^\ 3z 2 -r 2 

— In — + y/1 + (r r b ) 2 + c x ^+ 



x 2 — y 2 



rl(l + c 4 (r/r b ) 2 f 2 



(A2) 
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where r b and M b are the bulge core radius and the bulge mass respectively. Equation (A2) 
can be rewritten as a product of a dimensional factor and a dimensionless function as: 



U b (x',y',z') 



47T- 



GM h 



n 



— In [r' + y/l + r 



a 



II „/2 



3z' z - r 



<"i 



H- 



2(1 + c 2 r /2 ) 3 / 2 



-3c H 



x' 2 - y' 2 
(1 + c 4 r' 2 ) 3 / 2 



where 



.r 



y 



r = — , x =—, y = — , z = — . 
n n r b r b 



(A3) 



(A4) 



Also the first term on the right side of Eq.(Al) can be written as the product of a dimensional 
factor and a dimensionless one: 



N 



N 






(A5) 



with m'i = rrii/M b . 

Finally, once defined a dimensionless time 



t' 



(A6) 



being t c 



l/GMb) 1 ' 2 the bulge crossing time, Eq.Al may be written as: 



A? 



2^ 77T3 ( r * - r j) + 47rVf/ • 



dt' 2 



(A7) 



This implies that, once assigned the initial conditions r^(0), v^(0), the existence of a unique 
solution for the Eq. (A7) ensures that all the results obtained can be scaled in terms of the 
ratios r/r b , mjM b and t/t cross . 



B. The construction of mult i mass King model 



As described in King (1966), in a single-mass isotropic King model the phase-space 
stellar distribution function is given by: 



f(r,v) 



a 



exp 



E 



C 
- I - exp I — 

mcr z I \a z 



, if E < -mC 



(Bl) 



otherwise 
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where 



E 



-mv 2 + mip(r) 



(B2) 



is the energy of a star, ip(r) is the mean gravitational potential generated by the cluster and 
a is a normalization constant. The 'global' parameter C is related to the tidal radius r t by 
the implicit relation 

1>(r t )+C = 0. (B3) 

The mass density can be found integrating the distribution function f(r, v) over the velocity, 
obtaining an explicit relation for p as a function of the potential i/j: 



P 



f(r,v)4irv 2 dv 

I E<-mC 

A7imae c ^(2a 2 ?/ 2 



(B4) 



ip + C 



a" 



1/2 



+ 



+ 



TX 



ij + c 

47imae c l°\2a 2 yi 2 



erf 



1 



4> + C 



o A 



ip + C 



<7 Z 



3/2' 



(B5) 



-^Z7 + 



= kp(U) 



3/2 



(B6) 

(B7) 



where U = (ip + C)/a 2 is the dimensionless potential, k = Airmae ^ (2a 2 ) 3 / 2 and p is the 
dimensionless density, which explicitely depends only on U. Once assigned initial conditions 
for the potential ip and its derivative ip' in r = 0, the Poisson equation 

d 2 ^/dr 2 = 4ttG P , 

^'(0) = 

can be rewritten in terms of the dimensionless potential, in the form: 

d 2 U/df 2 = 9p(U)/p(0) = 9p(U)/p(0) 
U(0) = U 
U'(0) = 



where p(0) = p(Uo) and r = r/r core being 



9a 2 
4iiGp 



:bi 
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the King radius. 

Once assigned as initial parameters Uq, U' q , the Poisson equation can be integrated, ob- 
taining the dimensionless potential U(r), the dimensionless mass density p(r) and the tidal 
radius r t , being r t = r t /r core with r core yet not determined. To determine the core radius 
r core and the costant k (which depends, among others, on the normalization constant a), it 
is possible to procede as follows. 

Once assigned as initial parameters the total mass of the cluster M tot and the velocity dis- 
persion a of stars in the system, using the following relations 



M, 



lot 



4np(r)r 2 dr 



Ankrl 



p{f)r 2 dr 



(B9) 



and the Eq.B8, it is possible to calculate r core and k and hence to obtain p(r), r core and 

In a multi-mass isotropic King model, as described in Da Costa & Freeman (1976), 
stars are first grouped in n different mass classes, each characterized by a mass Wj. The 
phase-space stellar distribution function for the i-th mass class is given by: 



fi(r,v) 



(•>•; 



cxp 



E; 



micrj 



C 
0X1,1 a 2 



, if Ei < —rriiC 



(BIO) 



otherwise 



where 



Ei 



-rriiV 2 +rriiip(r) 



(bh; 



is the energy of a star in the i-th mass class, ip(r) is the mean gravitational potential generated 
by the whole cluster, at is a normalization constant and C is related to the cluster tidal radius 
r t by Eq.B3. Once again, the mass density for the i-th mass class can be found integrating 
the distribution function fi(r,v) over velocities, obtaining an explicit relation for pi in terms 
of the dimensionless potential U defined above and the ratio a 2 /a 2 : 



Pi 



fi(r,v)4:7iv 2 dv 

Ei<—rriiC 

A7im l a l e c ^(2a 2 f 2 



(B12) 



4> + C" 



1/2 



07 



+ 



cxp 



ip + C 



o~t 



erf 



a 2 



+ 

1 
3 



±+c_y /2 ' 
' -i ) 



(B13) 
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^m l a l e c l^{2a 2 l fl 2 -- 



+ 



.a 



exp —U— I erf 



a; 




-U 



a 
a 2 



2\ 3/2- 



= k t p t (U,a 2 /a 2 ) 



(B14) 
(B15) 



where ki = Anm i aie c / (T ^ (2a 2 ) 3 / 2 and pi is the dimensionless density. Note that the density 
profiles pi of the % — th mass class are related to the "global" density distribution p by the 
relation: 



p( r ) = ^2pi( r )- 



(B16) 



i=i 



To distribute stars in the cluster according to this isotropic multimass King model, we 
proceeded in the following way: 

• Once assigned U , U' Q , the total cluster mass M tot and the velocity dispersion <x, we 
integrate the Poisson equation as in the case of a sigle mass model previously described, 
obtaining the dimensionless potential U(f), the "global" mass density p, the core radius 
r core and the tidal radius r t . 

• Then we assigned the mass rrti of stars in the i-th mass class and the total mass M to t,i 
of each mass class (i.e. M totji = n; x toj, being n» the number of stars populating 
the i-th mass class). In our case, we chose to set the masses M to t,i according to the 
Salpeter's mass function. For a given value of the ratio a 2 jo\ (once obtained all the 
other values according to energy equipartition using the relation mi<j 2 = mi<j 2 for 
i > 2), we calculate pi(U,a 2 /a 2 ) and then the coefficients ki using Eq.B9 applied to 
the i — th mass class. 

• We varied the ratio a 2 /a 2 until the relation B16 is satisfied with the desired accuracy. 
Finally, stars velocities were generated according to Eq.BlO. 



C. The GC initial mass function 



The GC we considered in our simulations is supposed to have an age of t 



<)<■ 



10 



0.5 



yr, thus only stars more massive than ~ 1.2 M Q are at that time evolved up to a compact 
remnant (Straniero et al 1997; Dominguez et al 1999). For this reason, we considered only 
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stars distributed according to the Salpeter's MF with masses in the range 0.12 < m < 1.2 
M Q . Moreover, we assumed that the contribution to the low mass population due to the 
mentioned remnants is practically negligible. 

Indeed, according to the Salpeter's MF, the ratio between the number of remnants 
whose progenitor had a mass around m p and the number of stars with mass around m is 



N / x 2 - 35 



N m \ m p 



(CI) 



Supposing that such progenitors are those giving rise to remnants with mass m, then, from 
the estimates in Straniero et al (1997) and Dominguez et al (1999), m p (M Q ) ~ 9.5(m— 0.45), 
with m > 0.45 M Q because stars with lower mass remnants cannot be evolved in a Hubble 
time. Thus, substituting in Eq. (CI), 



-/Vremn /0.11 X 171^ 



rs_/ 



N m \m- 0.45 

~ otherwise 



if m>mi, (C2) 



where mi is the lowest mass a remnant can have at the assumed cluster age. From fitting the 
above-cited estimates, this lower limit turns out to be mi ~ (0.451ogt 9C — 1.2)~ 31 + 0.45 ~ 
0.59 M . 

One can see that the ratio in Eq. (C2) is monotonically decreasing for m > mi, hence 
the maximum takes place for the lowest mass class we used in the model, i.e. m — 0.71 M , 
for which iV remn /JVo.7i ~ 0.05. Since in our numerical representation N ji/N ~ 0.01 (N is 
the total number of particles), then in this class there should have been about 5 x 10~ 4 iV 
remnants. Thus, bearing in mind that the less populated mass class contains ~ 2 x 10~ 3 iV, we 
can reasonably affirm that the MF we assumed for the initial conditions was not substantially 
affected by stellar evolution neither at the initial time t gc nor later during the simulation 
(because it lasts much shorter than t gc ). 
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Fig. 1. — First orbital period of a 3 x 1O 5 M globular cluster in the potential described in 
Sec. 2.2. The cluster moves on a quasi-circular orbit (e = 0.03) around the galaxy center in 
a clockwise direction (see text). Distances are in units of the galactic bulge radius r&. Some 
snapshots are labelled with time, expressed in units of the galactic bulge crossing time t cross . 
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Fig. 2. — Second orbital period of the globular cluster described in the Fig.l. 
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Fig. 3. — Third orbital period of the globular cluster described in the Fig.l. 
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Fig. 4. — Last orbital period of the globular cluster described in the Fig.l. 
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Fig. 5. — First orbital period of the 3 x 1O 5 M globular cluster in a loop orbit with ellipticity 
e = 0.27 around a triaxial galaxy. The cluster moves in a clockwise direction. Some snapshots 
are labelled with time. The dotted line represents the cluster orbit. 
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Fig. 6. — Second orbital period of the globular cluster described in the Fig.5. 
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Fig. 7. — Third orbital period of the globular cluster in the Fig.5. 



-28- 



i o 

N 



i 1 r 



~i 1 1 r 



t = 28.8 t 



J I I L 



t 1 1 r 



t = ... 



# 



.''2W" 



J I I L_ 



-5 





y[ r J 



i 1 r 



Fig. 8. — Snapshots of the 3 x 10 5 M Q globular cluster in a loop orbit with ellipticity e = 0.57 
around a triaxial galaxy. The cluster moves in a clockwise direction. Some snapshots are 
labelled with time. The dotted line represents the cluster orbit. 
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Fig. 9. — Snapshot of the 3 x 1O 5 M globular cluster in the loop orbit with e = 0.27 at 
t=23.1t cross . The upper panel shows the system and the orbit described by the cluster density 
center (solid line). It is evident from the zoom in the bottom panel that the tails around the 
cluster core can lead to not reliable information about the orbital path of the cluster (solid 
line). 
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Fig. 10. — Volume mass density of the cluster in the case of the quasi-circular orbit (e = 0.03), 
at four different epochs, as labelled. The dashed line in each panel represents the best King 
model fit at that epoch. The presence of clumps in the tails are clearly visible. 
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Fig. 11. — Volume mass density of the cluster moving on the loop orbits with e = 0.27 (left 
column) and e = 0.57 (right column), at different epochs. The dashed line in each panel 
represents the best King model fit at that epoch. Note that in the case of most eccentric 
orbit, clumps are not yet formed at t=19.2t cross . 
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Fig. 12. — Surface density profile of the cluster at t=28.8t cross in the case of the quasi-circular 
orbit. Two different regions are plotted: that containing the trailing tail (filled triangles) 
and that containing the leading tail (open squares). The line-of-sight is parallel to the x axis 
and so perpendicular to the orbital plane. 
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Fig. 13. — Panel (a): Linear mass density as a function of the curvilinear abscissa s, set 
equal to zero at the cluster center, negative for the traling tail and positive for the leading 
tail. Panel (b): Derivative of the stellar tangential velocity respect to s. 
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Fig. 14. — Stars velocities along the x axis vs. the distance from the cluster center, for the 
case of the quasi-circular orbit (e = 0.03), plotted at four different epochs. Once stars begin 
to escape from the cluster, the velocity profile shows a minimum and then increases again. 
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Fig. 15. — Velocity dispersion profiles along the x axis for the cluster in the quasi-circular, 
at four different epochs. 
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Fig. 16. — Stellar mass belonging to the cluster, in units of the initial cluster mass, as a 
function of time, expressed in units of the bulge crossing time. Solid line: quasi-circular 
orbit with N — 1.6 X 10 5 particles. Dashed line: quasi-circular orbit with N — 1.6 x 10 4 
particles. Dotted line: loop orbit with e = 0.27. Dot-dashed: loop orbit with e = 0.57. 
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Fig. 17. — Left column: time evolution of the mean mass of stars in three different regions of 
space centered with the cluster. Panel (a): Mean stellar mass inside r = 0.016rb, in the case 
of quasi-circular orbit (solid line), loop orbit with e = 0.27 (dashed line) and loop orbit with 
e = 0.57 (dot-dashed). Panel (c): Mean stellar mass between r = 0.016?";, and r = 0.036rb. 
Panel (e): Mean stellar mass outside r = 0. 0367V Right column: evolution of the mean mass 
of stars in three different regions of space as a function of the fraction of mass lost from the 
system (in this case both the mass lost and the mean stellar mass have been averaged on 
time interval of 2.9 t cross for the two loop orbits with greater ellipticities). Panel (b): Mean 
stellar mass inside r = 0.016r{„ in the case of quasi-circular orbit (solid line), loop orbit with 
e = 0.27 (open circles) and loop orbit with e = 0.57 (solid circles). Panel (d): Mean stellar 
mass between r = 0.016?";, and r = 0.036r{,. Panel (f): Mean stellar mass outside r = 0.036r&. 
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Fig. 18. — Mass function of stars belonging to the cluster at three different epochs: when the 
cluster has lost the 20% (solid line), the 35% (dashed line) and the 75% of its initial cluster 
mass (dot-dashed line). Only the evolution of the mass function in the case of quasi-circular 
orbit is shown, because the other curves (corresponding to orbits with greater ellipticities) 
coincide with these. 
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Table 1. Orbital parameters. 



N 




e 


(^0, 2/0: 


• z ) 


(V x ,0, 


v y ,o,v Z: o) 


1.6 x 


10 5 


0.03 


(0,0,2 


.50) 


(0, 


,1.82,0) 


1.6 x 


10 4 


0.03 


(0,0,2 


.50) 


(0. 


,1.82,0) 


1.6 x 


10 4 


0.27 


(0,0,3 


.50) 


(0. 


,1.30,0) 


1.6 x 


10 4 


0.57 


(0,0,7 


.50) 


(0: 


,0.78,0) 



Note. - - Orbital parameters of the cluster in 
the four simulations performed. The first col- 
umn shows the total number of particles used in 
each simulation. Initial positions and velocities 
of the cluster with respect to the galaxy center 
(columns 3 and 4) have been expressed, respec- 
tively, in units of the galaxy bulge radius r b and 
of the bulge typical velocity dispersion r b /t CTOSS 
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Table 2. Clumps emersion from the background. 



r[r b ] 


Pcl/P* 


''[''cross J 


2.04 


<0.1 


9.62 


3.44 


0.3 


13.5 


3.47 


0.3 


19.2 


2.12 


<0.1 


23.1 


3.44 


0.3 


25.0 


2.05 


<0.1 


28.8 



Note. - - Clumps emersion from 
the background density, in the 
case of loop orbit with ellipticity 
e = 0.27. In the first column the 
cluster distance from the galaxy 
center is given; the second col- 
umn shows the ratio between the 
clumps local density and that of 
the stellar background; the third 
column shows the time from the 
beginning of the simulation, ex- 
pressed in units of the bulge cross- 
ing time. 



